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BOUNDARY  VALUE  PROBLEM 
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n 


by 

L.  J.l  Richards 


SUMMARY 


This  Memorandum  describes  a Fortran  subroutine,  BOVA,  which  provides  the 
solution  of  a set  of  nonlinear  ordinary  differential  equations  subject  to  non- 
mixed  two-point  boundary  value  conditions.  A finite  difference  method  is  used, 
together  with  Newton's  method.  A specification  and  listing  of  the  subroutine  is 
given . 
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1 INTRODUCTION 

We  consider  the  following  problem.  Given  a set  of  n nonlinear  first 
order  ordinary  differential  equations 


Xj  = *2’  *n^ 

^2  ” ^ 2 1 * ^2*  • • * * ^n^ 

X — f (X.  } • • • > 

n n I z n 


(1-1) 


find  a solution  trajectory  x,(t),  x„(t) x (t)  over  the  interval  (T.,T  ] 

1 / n Ur 

which  satisfies  n two-point  boundary  conditions  each  of  the  form 


or 


"i«0> 


d. 

1 


"iWp’ 


d. 

1 


for  some  i(l  ^ i s n)  . The  differential  equations  are  autonomous,  ie  the 
independent  variable  t does  not  appear  in  the  right  hand  sides  of  (1-1).  A 
non-autonomous  system  can  easily  be  put  into  autonomous  form  by  defining  a new 
variable  x , (t)  = t , satisfying  x , = 0 , with  x , (T.)  = T„  . Similarly, 
equations  involving  higher  order  derivatives  can  be  replaced  by  a larger  set  of 
first  order  equations  by  defining  new  variables.  We  must  assume  however  that 
each  of  the  boundary  conditions  holds  at  only  one  point  (either  T^  or  Tp) . 

If  this  assumption  were  not  satisfied  a method  similar  to  that  described  in 
section  2 could  be  used,  but  this  would  involve  a large  increase  in  the  amount 
of  computation  and  storage  required. 

The  Fortran  subroutine  described  in  this  Memorandum  was  designed  for  use  in 
a numerical  application  of  optimal  control  theory.  When  an  optimal  control 
problem  is  solved  using  Pontryagin's  Maximum  Principle,  this  leads  to  a two-point 
boundary  value  problem.  It  was  found  that  this  boundary  value  problem  could  not 
be  solved  by  conventional  'shooting'  methods,  so  that  an  alternative  method  was 
required.  The  method  described  here  is  applicable  to  any  two-point  boundary  value 
problem,  but  it  would  not  be  a good  choice  for  a general  purpose  two-point 
boundary  value  algorithm,  since  in  most  cases  more  effective  algorithms  can  be 
found . 
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2 METHOD 

Equations  (1-1)  and  (1-2)  may  be  rewritten  in  a more  compact  notation. 

The  differential  equations  are 

X = f(x)  (2-1) 

and  the  boundary  conditions  are 


C°x(Tg)  = d° 
C^x(Tp)  = d^ 


(2-2) 


0 F 

where  x is  an  n-vector,  d is  an  H-vector,  d is  an  (n  - i.)-vector,  and 
0 ~ F 

C and  C are  i x n and  (n  - H)  x n matrices  with  constant  elements. 

We  use  a finite  difference  formula  to  convert  the  set  of  nonlinear  differential 
equations  (2-1)  into  a set  of  nonlinear  algebraic  equations,  which  are  then 
solved  using  Newton's  method.  Instead  of  trying  to  determine  the  solution 
trajectory  throughout  the  interval  we  consider  its  value  only  at  regularly  spaced 
'grid  points' 


T = T„  + rh 
r 0 


(r  = 0,...,N) 


where  h = (T  - T )/N  and  N is  some  positive  integer.  The  method  will  yield 
F U 


x(Tq),  x(T|) 


The  finite  difference  approximation  is  similar  to  one  suggested  by  Fox 


y,  - Vq  = ih(y,  + y^)  + Cy^  . 


However,  we  omit  the  correction  operator  Cy^  and  use  the  trapezoidal  rule 


Fl  - Fq  = ih(yj  + y^)  . 


(2-3) 


This  leads,  in  our  notation,  to  the  vector  difference  formula 


x(T^^|)  - x(T 


^)  = ih(f(x(T^^,))  + f(x(T^))j 


(r  -0,..  ..N-1) 
(2-4) 
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It  is  convenient  to  replace  the  N + 1 vectors  x(T^)  each  having  n 
components  by  a single  vector  z having  n(N  + 1)  components.  Let 


where 


k = i + rn 

(i-  l,...,n,  r-0,...,N).  (2-5) 


The  equations  (2-4)  provide  nN  equations  involving  the  components  of  z . 

A further  n equations  come  from  the  boundary  conditions  (2-2),  the  first  I 

of  these  equations  involving  the  components  Zj z^  and  the  remaining 

(n  - 1)  equations  involving  the  components  » •••>  ^^N+n  ’ thus  have  a 

total  of  n(N  + 1)  nonlinear  algebraic  equations  in  n(N  + 1)  unknowns,  and  we 
can  write  these  equations  as 


g(z)  = 0 


(2-6) 


where  8j^(z) 


^ ‘^k2"2 


+ c,  z - d, 
kn  n k 


(k  = 1 H) 


(2-7) 


8k 


^rn+i  ^(r-l)n+i 


- ih(f.(z  ^ z . ) + f.(z.  , ,,  z,  ^ )) 

\ 1 rn+i  ’ rn+n  i (r-l)n+l  (r-l)n+n/ 

(k=i+rn+i.,  i = l r *=  0, . . . ,N-1 ) (2-8) 


8k  (z) 


F F F F 

c • I 2 , + c . „z  +'...+  c . z - d . 

nN+l  j2  nN+2  jn  nN+n  j 

(k»nN+t+j,  j = 1 , . . . ,n-S.) . (2-9) 


2 3 

Newton's  method  for  solving  nonlinear  equations  is  well-known  ’ and  will 

be  only  briefly  described  here.  Starting  with  an  initial  guess  z^  the  method 

I 2 

generates  successive  approximations  z , z , ...,  to  the  true  zero  z . The 
iterative  process  is  defined  by 


i+l 


z^  - J“'(z^)(g(z^)) 


(2-10) 
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where  J(z^)  is  the  Jacobian  matrix  evaluated  at  , ie  the  matrix  whose 

(p,q)th  element  is  (3g  /9z  ) (z^)  . If  ^ is  sufficiently  close  to  z the 

0 12^*^  ” . 

sequence  5 » z , z , converges  quadratically  to  z . 

3 COMPUTATIONAL  DETAILS 

In  using  equation  (2-10)  it  is  not  necessary  to  find  the  inverse  of  the 
Jacobian  explicitly.  It  is  sufficient  to  determine  the  vector  satisfying 

J(z^)e^  = g(z^) 

and  then  put 

i+1  i . 

z “ z + 

The  solution  of  the  set  of  linear  equations  (3-1)  is  made  easier  by  the  fact  that 
the  matrix  J(z^)  has  most  of  its  elements  equal  to  zero.  In  fact  it  has  the 
form  shown  in  Fig  1 . This  means  that  use  can  be  made  of  an  algorithm  which  takes 
advantage  of  the  band  form  of  the  matrix  to  reduce  the  amount  of  storage  and 
computation  time  required. 

In  order  to  calculate  J(z)  one  must  be  able  to  calculate  partial 

derivatives  of  the  functions  f,,  f_.  ...,  f with  respect  to  each  of  -he 

I 2 n 

variables  Xj , x^,  . . . , x^  . It  may  be  that  these  partial  rivatives  are 
difficult  or  impossible  to  calculate  directly,  for  example  the  functions  ai 
given  not  in  analytical  form  but  as  a series  of  tabular  values.  In  this  case 
it  will  be  necessary  to  make  an  approximation  to  the  derivatives  by  some 
numerical  means.  In  Appendix  B a subroutine  is  listed  which  does  this  using  the 
simplest  of  finite  difference  approximations. 

The  success  of  Newton's  method  depends  on  the  distance  of  the  initial 
guess  z*^  from  the  true  solution  z . The  larger  this  distance  is,  the  greater 
the  number  of  iterations  needed  to  converge  to  a given  accuracy.  If  the  distance 
is  too  large  the  algorithm  may  not  converge  at  all,  and  the  method  will  'blow-up', 
ie  produce  estimates  which  increase  in  magnitude  rapidly  and  without  limit. 

4 SPECIFICATION  OF  BOVA 

Title  Solution  of  two-point  boundary  value  problem  using  a finite 
difference  method. 


(3-1) 

(3-2) 


Programing  language  1900  Fortran. 
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Description  The  subroutine  produces  an  approximate  solution  to  the 
differential  equations  (1-1)  subject  to  boundary  conditions  (1-2).  The  solution 
is  given  at  a finite  number  of  regularly  spaced  grid  points. 

Accuracy  Single  precision  arithmetic  is  used,  and  the  accuracy  obtained 
depends  on  the  number  of  grid  points  employed  and  on  the  input  parameter  ACC. 

Experience  The  subroutine  has  performed  satisfactorily  on  a number  of 
test  cases.  See  section  5. 

Method  See  sections  2 and  3. 

The  call  statement  The  CALL  statement  is 

CALL  B0VA(M1,  M2,  M3,  M4,  M5,  M6,  M7,  MAXIT, 

HH,  ACC,  ICND,  Z,  NIT,  EPS,  AA,  AM, 

INT,  XY,  F,  FY,  P) 


where  Ml,  M2,  M3,  M4,  M5,  M6 , M7,  MAXIT,  NIT  are  integer  variables. 

HH,  ACC  are  real  variables.  ICND  is  an  integer  array  of  length  at  least  Ml. 

INT  is  an  integer  array  of  length  at  least  M3.  AA  is  a real  array  of  length  at 
least  M4.  AM  is  a real  array  of  length  at  least  M6 . XY,  F,  FY,  P are  real  arrays 
each  of  length  at  least  Mi  , 

Input  parameters  Before  calling  the  subroutine,  the  user  must  set  the 
arguments 

Ml  = n,  the  number  of  equations 

M2  = N + 1,  where  N determines  the  number  of  grid  points 


M3  = Ml  X M2 
M4  = M3  X (3n  - 1 ) 


M5  = Jl,  the  number  of  boundary  conditions  specified  at  the  initial 
point  Tq 

M6  - M3  X (m5  + Ml  - 1) 

M7  = reporting  parameter.  If  M7  > 0 the  subroutine  REP  is  called  after 

every  M7  iterations.  If  M7  ^ 0 REP  is  not  called.  See  'subroutines 
called ' . 


MAXIT  = the  maximum  number  of  iterations  allowed 

HH  = h,  the  distance  between  grid  points  = (T„  - T„)/N 

r U 

ACC  = accuracy  required.  The  iterative  process  is  stopped  if  the 
norm  of  the  vector  difference  between  successive  iterates  is 
less  than  ACC.  The  norm  used  is  the  uniform  or  Chebyshev  norm: 


c||  e^l : i = 1 , . ...  m| 


lCND(i)  = (i  = 1 , . . .,  Ji) 

= (i  = Jl  + 1 , . . . , n) 

where  Xj^  is  the  ith  variable  specified  at  the  initial 
point  Tq  , and  xj^^  is  the  (i  - Jl)th  variable  specified  at 
the  final  point  Tp,  . 

Results  At  exit  we  have 

NIT  = the  number  of  iterations  needed  for  convergence  to  the  required 

accuracy.  If  MAXIT  iterations  were  performed  without  convergence 
to  the  required  accuracy  NIT  is  set  to  -I . If  the  iterative 
process  'blows  up',  ie  produces  estimates  which  increase  in 
magnitude  rapidly  and  without  limit  NIT  is  set  to  -2. 

Z(k)  = x^(T^)  (k  = i+rn,  i = l,...,n,  r = 0,...,N). 

State  of  other  variables  All  the  input  parameters  are  left  unchanged. 

The  arrays  EPS,  AA,  AM,  XY,  F,  FY,  P are  used  as  work  areas. 

Subroutines  called  The  user  must  provide  the  following  four  subroutines: 

(1)  Subroutine  GUESS (Z ,M,HH) 

On  entry,  this  subroutine  must  accept  the  values  of  M = n(N  + 1)  and 
UH  = (Tp  - ^0^^^  termination  must  provide: 

Z(k)  = x?(T^)  (k  = i+rn,  i = I r = 0,...,N] 

where  ^ guessed  approximation  to  x^(T^)  . Z is  a real  one- 

dimensional array  of  length  M . 


(2)  Subroutine  FUN(X,F,M) 

On  entry  this  routine  must  accept  the  values  of  X = x and  M * n , and 
on  termination  must  provide: 
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F(i)  = f^(x)  (i  - 1 , . . . ,n) 

X and  F are  real  one-dimensional  arrays  of  length  M . 

(3)  Subroutine  PARD(X, I,P,M) 

On  entry  this  subroutine  must  accept  the  values  of  X = x,  I = i and 
M = n , and  on  termination  must  provide 

3f. 

P(j)  “ (x)  (j  = 1 n)  . 

j 

X and  P are  real  one -dimensional  arrays  of  length  M . 

(4)  Subroutine  REP(Z,M,AN0RM,IR) 

This  subroutine  enables  the  user  to  monitor  the  progress  of  the  iterative 
process.  On  entry 

M = n(N  + 1) 

IR 

Z(k)  * , the  IRth  approximation  to  x^(T^) 

(k  “ i+rn,  i = l,...,n,  r = 0,...,N) 

ANORM  = = ||Z^^  - 

where  ||e|l^  = max|le^l:  i = 1,  m| 

Z is  a real  one -dimensional  array  of  length  M . 

The  ICL  subroutines  FPDEBM  and  FPMOVE  are  also  used. 

5 TEST  RESULTS 

We  consider  the  problem 

••3  2 

y • y - sin  t(l  + sin  t) 

over  the  interval  0 to  r , with  boundary  conditions 

y(0)  ” y(n)  » 0 


which  has  the  solution 
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Expressing  the  problem  in  the  form  of  equation  (1-1)  it  becomes 


3 2 

Xj  - sin  x^(l  + sin  x^) 


Xj(0)  = 0 


x^Ctt)  = 0 


X2(0)  = 0 


J 


The  numerical  results  obtained  will  depend  on  the  number  of  grid  points, 
the  accuracy  demanded  and  the  initial  guess.  The  following  table  gives  the 
results  of  a number  of  runs. 


Case  No . 


Maximum  error 


No.  of  iterations 


1 

2 

3 

4 

5 

6 


0.0003 

0.002 

0.00005 

0.0003 

0.0003 

0.0003 


5 

5 

5 

6 
6 
9 


Case  1 had  N = 40,  ACC  = 0.001  and 


x"(T_^) 


0 

0 

T 


> r = 0,  1 N 


Case 

2 - as 

for 

case 

1 . 

except 

N = 

20 

Case 

3 - as 

for 

case 

1, 

except 

N = 

100 

Case 

4 - as 

for 

case 

1, 

except 

ACC 

= 0.00001 

Case 

5 - as 

for 

case 

1. 

except 

the 

subroutine  PARD  listed 

was 

used  to 

approximate  the  partial  derivatives 

Case 

6 - as 

for 

case 

1 , 

except 

the 

initial  guesses  were: 

A 


with  boundary  conditions 


Xj(0)  = x^(0)  = x^(0)  = 1 

x^(0)  = 0 

Xj(n/^)  = x^(ti/4)  = X|q(it/4)  - - /z/2 

x^(it/4)  - - I 

x^(ir/4)  - 1 


Xg(ir/4) 


0 
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The  exact  solution  is 

x^(t)  = cos  it  (i  “ 1,...,5) 

x^(t)  = sin  (i  - 5)t  (i  = 6,. ..,10). 

The  equations  were  solved,  using  subroutine  BOVA,  with  N = 40,  ACC  = 10  and 
initial  guesses  x?(T^)  =0.5  for  all  values  of  i and  r except  those 
specified  in  the  boundary  conditions.  Six  iterations  were  needed  for  convergence. 
The  maximum  error  was  0.001,  this  error  occurring  in  Xg  , which  was  not 
specified  by  the  boundary  conditions  at  either  end  of  the  interval.  On  the 
ICL  1906S  computer  this  run  took  43  seconds. 


i 
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SUBmtJTIME  BOVACM  1 ACC 

l,ICMD,Z,niT>EPS,AA,AM>lNT,XY,F,FY,P) 

INTECER  INT(M3) 

REAL  Z(M3) ,EPS(M3) , AA(MA) / AM(M6)^XY(M1 ) /F<M! )/ FYCMl ) , P< 
IMl  ) 

INTEGER  ICNDCMl) 

COMMON/D1M5/N>NN^MNN/NDI,ND2,NE,  I STA>  1 P.P,H,H2 

N=M1 

NN  = M2 

NNN=M3 

NDI=Mi:j 

ND2=M6 

1STA=M5 

1RP=M7 

NE  = M']  '1-N 

H=HH 

H2=0.5»H 

CALL  BVPACMAXIT, ACC,NIT,Z,EPS, ICND, AA, AM,XY,F,FY 
UP,  IMT) 

RETURN 

END 


i J 


no  o n o<ooo 


41 
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Appendix  A 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

nvIDED 


SUDEDIITINF.  BVPA<  MAX  I T,  ACC , M I T, 7- , EPS>  I CM D,  AA,  AM,XY,  F,  FY,  P,  INT) 

THIS  SUDP.nUTltJE  SOLVES  A MONLIMEAR  TVO-POIMT  BOUNDARY  VALUE 
PROBLEM  USING  FINITE  DIFFERENCE  FORMULAE  AND  A 
NENTON-RAPHSO'J  PROCEDURE 


MAX  IT  IS  THE  MAXIMUM  NUMBER  OF  ITERATIONS  ALLOVED 

THE  PROCESS  STOPS  VHEN  TVO  SUCCESSIVE  ITERATES  AGREE  TO  AN 
ACCURACY  OF  ACC^OR  IF  THE  SOLUTION  BLOWS  UP 
NIT  IS  THE  NUMBER  OF  ITERATIONS  NEEDED 
(-1  IF  NO  CONVERGENCE, -2  IF  SOLUTION  BLOWS  UP  ) 

A CALL  IS  MADE  TO  SUBROUTINE  REP  AFTER  EVERY  I RP  I TERATI ON S, PR 


THAT  IRP  IS  +VE 

C CONTAINS  INFORMATION  ABOUT  THE  BOUNDARY  CONDITIONS 
THE  OTHER  ARRAYS  IN  THE  SUBROUTINE  LIST  ARE  NEEDED  BY 


IMPRO 


1/4 


AT  EXIT  Z CONTAINS  THE  SOLUTION 
INTEGER  R, INTCMNN) 

REAL  Z(NMU),EPS(NNN) ,AA(ND1 ) , AM( ND2 ) ,XY( N ) , F( N ) , FY< N ) , PC N) 
INTEGER  ICNDCN) 

CDMMON/DIMS/N>NN,NNN,NDl ,ND2,NE, I STA, IRP,H,H2 

FIRST  WE  MUST  MAKE  AN  INITIAL  GUESS  AT  THE  SOLUTION 
CALL  GUESS(Z,NNN,H) 

ANOrc-l=0. 

DO  1/4  I = 1,NNN 
/iN  = ADS<Z(  I > ) 

I F( AN  - AMORM ) 1 4,  1 4, 0 
ANORM=AN 
CONTINUE 
BL0W=AN0R:^*  1 . ES 


H 


i 


DO  3 P.=  1,MAXIT 

CALL  IMPRGVE(Z,EPS. I CND, AA, AM,XY, F, FY, P, INT) 
NEW  ITERATE  =QLD  ITERATE  ♦ EPS 
INFINITY  NORM  = ANORM 
ANOrM=0. 

DO  4 i=i,nn:; 

ZC I j=Zi I ;»EPSC I ) 

AN=ABS(EPS< 1 ) ) 

I FC  AN- ANORM)  4,4,0 
ANOPJ-!=AN 

4 CONTINUE 
IFCIRP)  8,8,0 
1FCR/IRP*IRP-R)  8,0,8 
CALL  REP<Z,NNM,ANOP-M,P) 

8 CONTINUE 

IFC  ANORM-ACC)  S,0,0 
IFC  AM0R;-J-BL0W)  3,0,0 
NIT=-2 
GO  TO  9 
3 CONTINUE 
NIT=-1 
RETURN 

5 NIT-R  ^ 

9 return 

END  i 


i 


non  ooonooono  n non 


Appendix  A 


15 


SUBRGUTIME  i;<!PrDVEC Z > EPS / I C‘JD,  AA,  A:-I>XY,  ?,  FY,  P,  IMT) 

THIS  SUBHOUTINE  TAKES  A GUESS  Z AMD  CALCULATES  EPS  SO  THAT 
Z + EPS  IS  A BETTER  APPP.OX  I MAT  I ON  TO  THE  SOLUTION  THAN  Z 
A NEVTOM-RAPMSON  METHOD  IS  USED 
INTEGER  R,  IMT(NNN) 

REAL  ZCMMM),EPS(MMM>>AA(ND1), AM(ND2),XY(N),F(N),FY(N), P(N) 

INTEGER  ICND(M) 

C0MM0M/DIMS/M,N:J,NMN>ND1,ND2,NE,  I STA,  IRP,H,H2 

FIRST  ZEROISE  THE  ARRAY  AA  WHICH  WILL  CONTAIN  THE  MATRIX 
AA< 1 )=0. 

CALL  FMOVE< AA< 1 ),AA( 2),MD1 - 1 ) 


R DEMOTES 
I DENOTES 
NNN=M*MN. 

NE=NMM-N 

ND1=(M1+M2+ 1 )*MMN 
NUMBER  OF  NONZERO 
BAND  MATRIX  A.  TO 


THE  NUMBER  OF  THE  PIVOT  POINT  ( 1 
THE  NUMBER  OF  VARIABLE  (1  TO  N) 
THIS  IS  THE  ORDER  OF  THE  MATRIX 


TO  NN) 


/ ND2  =MI*NNN  /WHERE  Ml  AND  M2  ARE 
SUDDIAGONALS  AND  SUPERDI AGOMAL S OF 
TAKE  ADVANTAGE  OF  THE  BAND  FORM  WE 


THE  NONZERO  ELEMENTS  OF  A IN  THE  1 -DI  MEN  SI  OMAL  ARFA'; 


THE 
THE 
STORE 
AA/  SO 


THAT 


POINT 


C 

C 

C 

C 

C 

c 

c 

c 
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A(I/J)  = AA( 1 *MMN*(M1 -( I -J ) ) 

WE  HAVE  Mi=ISTA+N-l  M2=2N-1-1STA 

WHERE  ISTA  IS  THE  NO.  OF  BOUNDARY  CONDITIONS  AT  THE  INITIAL 
H IS  THE  STEP-LENGTH.  H2=0.5*H 

THE  RIGHT  HAND  SIDES  OF  THE  EQUATIONS  ARE  STORED  IN  EPS 
THIS  WILL  BE  OVERWRITTEN  LATER 

FIRST  DEAL  WITH  BOUNDARY  CONDITIONS 

THE  DIFFERENCE  EQUATION  PROVIDES  (NN-n*N  EQUATIONS 
THE  REMAINING  N EQUATIONS  COME  FROM  THE  BOUNDARY  CONDITIONS 
SPLIT  BETWEEN  INITIAL  AND  FINAL  POINTS 
M1=ISTA+N- I 
M2=N+N-1-ISTA 
DO  A K=1/1STA 
EPS(K)=0. 
l=ICNDCK) 

AA(  K+N NN* ( M 1 ♦ I -K ) ) = 1 . 

CONTINUE 
DO  lA  K=ISTA+1/N 
IW I=NE+K 
EPS(IW1)=0. 

I=ICMD<K) 

AA(  IW  1 +NNN*  ( M 1 ■•■NE+  I - I W 1 ) ) = 1 . 

CONTINUE 


.y 
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C MOV  FDr.  THE  DIFFEPEtJCE  EQUATION 

C EACH  nOV  OF  A CONSISTS  OF  2N  NONZERO  ELEMENTS 

C THE  2ND  N COEFFICIENTS  OF  RO V< I STA+ C R- I >N+ I ) ARE  ALMOST 

C THE  SAME  AS  THE  1ST  N OF  RON  (ISTA+RN+I) 

C THE  COEFFICIENTS  INVOLVE  PARTIAL  DERI VATI VE S, AND  THE  POINTS 

AT 

C WHICH  THESE  ARE  EVALUATED  ARE  INDEPENDENT  OF  I / SO  VE  CAN 

C FIND  THESE  POINTS  IN  THE  OUTER  LOOP 

C 

C INITIALISE  FY  BEFORE  STARTING  OUTER  LOOP 

DO  12  K=UN 

12  XY<K)=ZC<) 

CALL  FUNCXY,FY,N) 

C SET  UP  1ST  N COEFFICIENTS  IN  FIRST  N EQUATIONS 

31  DO  22  I = UN 

CALL  PARDCXY,  UP,N) 

IRA=ISTA+1 

IU=MI-1RA 

C WE  NOW  FIND  A(1RA,K)  (K=1,...,N) 

DO  23  K=1,N 

23  AA<  IRA+NflN=K<  IW  + X)  > = -H2*P(X) 

IW 1=IRA+NNN*( IW+I ) 

AA( IWl)=AA( IWl )-l. 

22  CONTINUE 

DO  1 R=l,NN-2 

M=R*N 

L=M-N 

DO  3 I=1,N 

XY( 1)=Z(M+I  ) 

C ALSO  NOW  FIND  FUNCTION  VALUES 

3 F<  I ) = FY(  1 ) 

CALL  FUN(XY,FY,N> 

C BEGIN  INNER  LOOP 

32  DO  2 1=1, N 
IP.A=ISTA+(R-  1 )*>]+! 

C FIND  PARTIAL  DERIVATIVES  FOP.  2ND  N COEFFICIENTS  OF 

C EQUATION  IP.A  AiND  1ST  N COEFFICIENTS  OF  EQUATION  I RA+N 

CALL  PARDCXY, I,P,N) 

C WE  NOW  FIND  A(IPA,M  + ;<)  AND  A(IRA+N,M  + X) 

1W=M1+M-IRA 
DO  7 K=1,N 
IW1=IFA+NNN*( IW+K) 

1W2=IW  UN+(  1-NNN) 

AA( IWl )=-H2*P(K) 

7 AA( IW2)=AA( IW 1 ) 

IW1=IRA+NNN*< IW+I ) 

IW2=r.»UN*(  1-NMN) 

AAC IW I )=AA< IWl )* 1 . 

AA< IW2)=AA< 1W2)-1. 


i. 


i. 


} 
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C NOV  FIND  EPS(  I BA) 

EPS( 1BA)=H2*( F< I )+FY( I ) )-Z (M+1 )+Z(L+I) 

2 CONTINUE 
1 CONTINUE 

C WE  HAVE  MISSED  OUT  THE  2ND  N COEFFICIENTS  IN  THE  LAST  N EQU 

ATIONS 

C SO  SET  THESE  UP  NOW 

M=NE 
L=M-N 

DO  18  1=1, N 
XY( I )=Z(M+I ) 

18  F<I)=FY<I) 

CALL  FUN(XY,FY,N) 

33  DO  19  1=1, N 
IBA=I STA+L+I 
CALL  PARDCXY, I,P,N) 

IW=M1*M-IRA 
DO  20  K=1,N 
IWl  = inA+NNNt'(  I’/  + K) 

20  AA(  IWl  >=-H2'ttP(K) 

I W 1 = I n A+N  NN* ( I W + I ) 

AA( IW 1 )=AA( IW 1 )+ 1 . 

EPS<  lRA)=H2t'(  F<  I )+FY(  I ) ) -Z  ( M+ 1 ) +Z  ( L+ I ) 

19  CONTINUE 

C WE  HAVE  FINISHED  SETTING  UP  THE  EQUATIONS  AND  CAN  NOW  PRGCE 

ED 

C TO  SOLVE  THEM 

3A  1W=1 

C AM  AND  INT  ARE  WORK  ARRAYS  USED  BY  FPDEBM 

CALL  FPDEBM (NNN, IW,M1 ,M2, AA( 1 ), AM( 1 ) ,EPS< 1 ) , INT( 1 ) , I V2) 

C EPS  HAS  BEEN  OVERWRITTEN,  Af>JD  NOW  CONTAINS  THE  REQUIRED  SOL 

UTION 

CHECK  FOR  SINGULARITY 

IF<r-72)  0,10,10 

WRITE(2,  1 1 ) 

11  FORMATC IX, ! 5H  SINGULAR  MATRIX) 

21  FORMAT ( IX, 6E 15. 6) 

STOP 

10  RETUP.N 
END 


If  O'-', 
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Appendix  B 

LISTING  OF  SUBROUTINE  PARD 


Sl  I.-nOLTUh 
P.LAL 

PLAL  f(IC) 

Cu.i.lO.J  /LST/'OC  10>  10) 
li-  ( I-  1 ) 5/0/5 
CALL  Ft.J(X/?/.J) 

LO  1 J=l/.'l 
LLL  = C.  0 !*>•(  J) 

Ir  (ALSCLLD- 1 . L-8)  0/0/2 
LtL= I . L-2 
XX=X< J) 

X ( J)=XX+LLL 
CALL  FL.J(X/F/.4) 

X ( J )=XX 
DO  1 K=1/.’J 

0 ( ;{ / J ) --  ( F < K ) - ? < K ) ) / DLL 
DO  6 <J=  1 / 14 
P(  J)  = 0<  I/vJ) 

ni-TVPu4 

E.4D 
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